% MLEs and Test Statistics for Returns Data .. note needs to be run twice .. once for Largest and once for Smallest Returns
clear all
small = 1.0e-10;
pinv_tol = 1.0e-05;
big = 1.0e+8;

% Add new directory to path to find functions
addpath('../Matlab_Functions/');

% -- File Directories   
outdir = 'out/';
figdir = 'fig/';
matdir = 'mat/';
cv_dir = 'CVS_Files/';

% Set options for optimizers
options_fminunc = optimset('Display','off');
options_fminsearch = optimset('Display','off');

% Identifiers for files 
% -------------------------- VW Returns -- largest --------------
application_str = 'Returns_Large';
data_fname = [matdir 'Returns_Large.mat'];
load(data_fname);
Y_data = ret;
% -------------------------------------------------------
% -------------------------- VW Returns -- smallest --------------
% application_str = 'Returns_Small';
% data_fname = [matdir 'Returns_Small.mat'];
% load(data_fname);
% Y_data = ret;
% -------------------------------------------------------
Y = Y_data;
T = size(Y,1);
k = size(Y,2);

% Starting Values to use for MLE 
xsi_start = 0.5;
sigma_start = std(Y(:));
mu_start = mean(Y(:));

% Find MLE from GEV distribution: Order of parameters is xsi, sigma, mu 
x_start = [xsi_start sigma_start mu_start]; % Starting Values
f_unc = @(x) minus_llf_gev(x,Y);
[x_1,fval_1]=fminunc(f_unc,x_start,options_fminunc); 
[x_2,fval_2] = fminsearch(f_unc,x_start,options_fminsearch); 
  if fval_1 < fval_2
      x_mle = x_1;
      fval = fval_1;
  else
      x_mle = x_2;
      fval = fval_2;
  end
llf_unconstrained = -fval;
xsi_mle = x_mle(1); 
sigma_mle = x_mle(2);
mu_mle = x_mle(3);

% Compute Score and Hessian
[score_mat_unc,info_op_unc,info_h_unc] = score_info_GEV_xi_sigma_mu(Y,xsi_mle,sigma_mle,mu_mle);
score_mat_unc = score_mat_unc - mean(score_mat_unc);    % score_mat_unc should have mean zero ... impose this

T = size(Y,1);
info = info_h_unc;
cov_mle = inv(info)/T;
se_mle = sqrt(diag(cov_mle));
se_xsi_mle = se_mle(1);
se_sigma_mle = se_mle(2);
se_mu_mle = se_mle(3);

% Compute Nyblom Statistics -- test whether all the parameters are constant
L_Stat = compute_nyblom(score_mat_unc,info_h_unc);
pv_L_Stat = pvalue_LStat(L_Stat,length(x_mle));  % Traditional p-value -- regular asymptotics
% Read in parameters for adjustment factor
fname = [cv_dir 'cvx_test0_' num2str(k) '_' num2str(T) '.csv'];
c_parm = readmatrix(fname);
adj = exp(c_parm(1) + c_parm(2)*xsi_mle + c_parm(3)*xsi_mle^2);
adj_L_Stat = adj*L_Stat;

% Construct Test that mu = sigma/xsi
x_start = [xsi_mle sigma_mle];
f_con_mu = @(x) minus_llf_gev_con_mu(x,Y);
[x_1,fval_1] = fminunc(f_con_mu,x_start,options_fminunc);
[x_2,fval_2] = fminsearch(f_con_mu,x_start,options_fminsearch); 
if fval_1 < fval_2
      x_mle_constained_mu = x_1;
      fval = fval_1;
  else
      x_mle_constained_mu = x_2;
      fval = fval_2;
end
llf_constrained_mu = -fval;
%  Compute adnjustment factor for test statistic
fname = [cv_dir 'cvx_test3_' num2str(k) '_' num2str(T) '.csv'];
c_parm = readmatrix(fname);
adj = exp(c_parm(1) + c_parm(2)*xsi_mle + c_parm(3)*xsi_mle^2);
llf_dif = llf_unconstrained-llf_constrained_mu;
Stat_con_mu = adj*llf_dif; 

% Construct Test that mu = sigma/xsi and xsi = 1
x_start = [sigma_mle];
f_con_xsi_mu = @(x) minus_llf_gev_con_xsi_mu(x,Y);
[x_1,fval_1] = fminunc(f_con_xsi_mu,x_start,options_fminunc);
[x_2,fval_2] = fminsearch(f_con_xsi_mu,x_start,options_fminsearch);
if fval_1 < fval_2
      x_mle_constained_xsi_mu = x_1;
      fval = fval_1;
  else
      x_mle_constained_xsi_mu = x_2;
      fval = fval_2;
 end
llf_constrained_xsi_mu = -fval;
%  Compute adnjustment factor for test statistic
fname = [cv_dir 'cvx_test4_' num2str(k) '_' num2str(T) '.csv'];
c_parm = readmatrix(fname);
adj = exp(c_parm(1) + c_parm(2)*xsi_mle + c_parm(3)*xsi_mle^2);
llf_dif = llf_unconstrained-llf_constrained_xsi_mu;
Stat_con_xsi_mu = adj*llf_dif; 

% Construct Confidence Interval for xsi;
% Form CI for xsi
% Step 1 -- read in grid of values of xsi and critical values and construct expanded grid
%   fname = [cv_dir 'cvx_test1_' str_cvs '_' num2str(T) '.csv'];
fname = [cv_dir 'cvs_test1_' num2str(k) '_' num2str(T) '.csv'];
tmp = readmatrix(fname);
xsi_grid_file = tmp(1,:)';
c_grid_file = tmp(2,:)';
% Compute Grid of values for xsi for evaluation
xsi_min = min(xsi_grid_file);
xsi_max = max(xsi_grid_file);
n_xsi = 100;
xsi_grid = linspace(xsi_min,xsi_max,n_xsi);
critical_value_xsi_grid = interp1(xsi_grid_file,c_grid_file,xsi_grid)';
xsi_grid = xsi_grid';
% Step 2 construct log-likelihood value for each value in Xsi Grid
llf_vec_con_xsi = NaN(n_xsi,1);
x_start = [sigma_start mu_start]; % Starting Values
for i = 1:n_xsi
    xi_con = xsi_grid(i);
    f_con_xsi = @(x) minus_llf_gev_con_xsi(x,Y,xi_con);
    [~,fval_1]=fminunc(f_con_xsi,x_start,options_fminunc);
    [~,fval_2]=fminsearch(f_con_xsi,x_start,options_fminsearch); 
    if fval_1 < fval_2
        fval = fval_1;
    else
        fval = fval_2;
    end
    llf_vec_con_xsi(i)=-fval;
end
llf_dif = llf_unconstrained - llf_vec_con_xsi;
ii = llf_dif < critical_value_xsi_grid;
tmp = xsi_grid(ii==1);
CI_xsi_lower = min(tmp);
CI_xsi_upper = max(tmp);

% Construct a confidence interval for x_90 -- GEV 90th quantile
  % Get MLE
  x_90_mle = gevinv(0.9,xsi_mle,sigma_mle,mu_mle);
  % Adjustment factor
  fname = [cv_dir 'cvx_test2_' num2str(k) '_' num2str(T) '.csv'];
  c_parm = readmatrix(fname);
  adj = exp(c_parm(1) + c_parm(2)*xsi_mle + c_parm(3)*xsi_mle^2);
  
  % Step 1: Input a grid of values 
  q_min = 3.0;
  q_max = 6.0;
  n_q = 500;
  q_grid = linspace(q_min,q_max,n_q)';
  % Step 2 construct log-likelihood value for each value in Xsi Grid
  llf_vec_con_q = NaN(n_q,1);
  % x_start = [sigma_start mu_start]; % Starting Values
  x_start = [xsi_mle sigma_mle];
  for i = 1:n_q
    x_90_con = q_grid(i);
    f_con_x_90 = @(x) minus_llf_gev_con_x_90(x,Y,x_90_con);
    [x_1,fval_1]=fminunc(f_con_x_90,x_start,options_fminunc);
    [x_2,fval_2]=fminsearch(f_con_x_90,x_start,options_fminsearch); 
     if fval_1 < fval_2
         fval = fval_1;
         x_1 = x_1;
     else
         fval = fval_2;
         x_1 = x_2;
     end
     llf_vec_con_q(i)=-fval;
     x_start = x_1;
  end
  llf_dif = llf_unconstrained - llf_vec_con_q;
  stat = adj*llf_dif;
  ii = stat < 1.0;
  tmp = q_grid(ii==1);
  CI_x90_lower = min(tmp);
  CI_x90_upper = max(tmp);
  % Make sure these are not at bounds 
  if CI_x90_lower == q_min
      CI_x90_lower = NaN;
  end
  if CI_x90_upper == q_max
      CI_x90_upper = NaN;
  end

% % Save Results to output file
outfname = [outdir application_str '_LikelihoodResults.txt'];
fid = fopen(outfname,'w');
fprintf(fid,'Data: %s\n',data_fname);
fprintf(fid,'T = %3i \n',size(Y_data,1));
fprintf(fid,'k = %3i \n',size(Y_data,2));
fprintf(fid,'MLEs with Hessian-based SEs (not correct, but nice to see) \n');
fprintf(fid,'xsi_mle = %6.4f (%6.4f) \n',[xsi_mle se_xsi_mle]);
fprintf(fid,'sigma_mle = %6.4f (%6.4f) \n',[sigma_mle se_sigma_mle]);
fprintf(fid,'mu_mle = %6.4f (%6.4f) \n',[mu_mle se_mu_mle]);
fprintf(fid,'MLEs 90th quantile of GEV: %6.2f \n',x_90_mle);
fprintf(fid,'95 percent CI for Xsi: %6.2f to %6.2f \n',[CI_xsi_lower CI_xsi_upper]);
fprintf(fid,'95 percent CI for 90th GEV quantile: %6.2f to %6.2f \n',[CI_x90_lower CI_x90_upper]);
fprintf(fid,'Some Test Statistics (5 percent CVs are 1.0) \n');
fprintf(fid,'   Test that mu = sigma/xsi: %6.2f \n',Stat_con_mu);
fprintf(fid,'   Test that mu = sigma/xsi and xsi = 1: %6.2f \n',Stat_con_xsi_mu);
fprintf(fid,'   Test for constant parameters (Nyblom)  %6.2f \n',adj_L_Stat);

fclose(fid);